function [bwOut ] = cutworm( bwImage, sktpAll, pieces, cutLocation)
%CUTWORM Return BW masks for specified sections of the worm
%   Input arguments
%
%    bwImage - a binary image of a worm
%
%   sktpAll - an ordered list of all skeleton points in the worm
%
%   pieces = the number of pieces the worm is to be cut into
%
%   cutLocation - the location of each cut specified as the index location
   %   of the pieceNo + 1 points defining equal sections of the worm


 %eg, n pieces have n+1 points, including endpoints 
 numPoints = pieces + 1 ;
 
 %Get the list of points and their relationship in the original skeleton
  [  ~, sktpIndices ] = createSktp(sktpAll, numPoints);
  
  bwOut = bwImage;
  %Iterate thru cut locations
  for i= 1: length(cutLocation)
  
      %Get nine points around each cut location
      ptList = sktpAll( sktpIndices( cutLocation(i))-4:sktpIndices(cutLocation(i))+ 4,:);
      
      %Get the number of points in the list (should be 9)
      ptNo = size(ptList,1);
      
      %Get the midpoint, should be 5
      midptLoc = ceil(ptNo/2);

      %Fit a curve to the nine points
      p  = polyfit( ptList(:,1),ptList(:,2),1 );

      % get the slope of the perpendicular line
      m = - 1/p(1);

      %Select the midpoint of the series of points as one point on the 
      %perpendicular line
      A = ptList(midptLoc,:);

      %Get the y-intercept 
      b = A(1,2) - m *A(1,1);
      
      %Find a point 25 pixels along the line in both directions
      phi = atan(m);
      B1(1,1) = A(1,1) - 25* cos(phi);
      B1(1,2) = m* B1(1,1) + b;
      B2(1,1) = A(1,1) + 25 * cos(phi);
      B2(1,2) = m* B2(1,1) + b;

      %Get a set of subscripts for the line segment between the two pts
      %Note that there will be many duplicates
      rpts = transpose( round(linspace(B1(1,1),B2(1,1),1000)) );   
      cpts = transpose( round(linspace(B1(1,2),B2(1,2),1000)) );   
      linePts = unique( [rpts, cpts], 'rows');
      
      %% Get a list of all perimeter points
      
      %Remove all interior pixels to get BW image of only the perimeter
      Iperim = bwmorph(bwImage, 'remove');
        
      %Get the row and column of the perimeter points
      [row, col] = find(Iperim);
      perimPts = [row,col];
      
      %% Find the intersecting points between line and perimeter
      
      %Get 
     
      for q = 1:2
         
          if q > 1
              SE = strel('square',2);
              Iperim = imdilate(Iperim, SE);
              [row, col] = find(Iperim);
              perimPts = [row,col];
          end
     
            intersecting = ismember(linePts, perimPts,'rows');
            [ rowIndex, ~] = find(intersecting);
            interPt = linePts( rowIndex,:);  

            interPtImage = false(size(bwImage)) ;
            index = sub2ind(size(bwImage), interPt(:,1), interPt(:,2));
            interPtImage(index) = 1;
            cc = bwconncomp(interPtImage);
        
      
          %Select one pixel index from each connected component to serve as
          %the intersection point. 
          interPtIdx = zeros(1, cc.NumObjects);
          for k = 1:cc.NumObjects
             pixelList  = cc.PixelIdxList{k};
             interPtIdx(k) = pixelList(1);
          end

          %Get the subcripts for the indices
          [interPtRow,interPtCol] = ind2sub(size(bwImage),interPtIdx);
          interPt = [interPtRow',interPtCol'];
          if size(interPt,1) > 1
              break;
          end
      end
      
        
      %% Find the intersecting points nearest the cut point
        w = repmat(A, size(interPt,1),1);
      distance = distl2( w , interPt );
      
      [ H, idx] = sort( distance, 'ascend');
      interPt = interPt( [idx(1),idx(2)],:);
      
      %% Cut the worm with a segment between the two points
        
      bwOut = drawsegment(bwOut,interPt(1,:),interPt(2,:), 0);

      %Remove diagnonal (8-connectivity between two sides of the break
      %line to make sure the section is not connected
      bwOut = bwmorph(bwOut,'hbreak');
      bwOut = ~ bwmorph(~bwOut,'diag');
      
  end
  
  
  
end


